Średnia długość życia wydaje się prostą do przewidzenia zmienną dla krajów, posiadających dane statystyczne. Klasyfikacja na poziomie 0.95 R2 jest dostatecznie wysoka, aby wykorzystywać dane statystyczne do zwiększania średniej życia w krajach – np. poprzez inwestowanie pięniędzy bądź czasu w czynniki, które najbardziej na to wpływają – są to między innymi poziom nauczania i ilość zgonów spowodowanych wirusem HIV.
Użyte bibliotek przedstawiam poniżej. Są to biblioteki do operacji na macierzach, statystyczne, do uczenia maszynowego oraz do przedstawiania wykresów.
usePackage <- function(p) {
if (!is.element(p, installed.packages()[,1]))
install.packages(p, dep = TRUE, verbose = FALSE)
require(p, character.only = TRUE, quietly = TRUE)
}
packageList <- c("dplyr", "matrixStats", "tidyr", "caret", "ggcorrplot", "plotly")
print(packageList)## [1] "dplyr" "matrixStats" "tidyr" "caret" "ggcorrplot"
## [6] "plotly"Instalacja bibliotek obdywa się za pomocą skryptu, z wyłączeniem logu.
sapply(packageList, usePackage)## dplyr matrixStats tidyr caret ggcorrplot plotly
## TRUE TRUE TRUE TRUE TRUE TRUEPowtarzalność obliczeń jest zapewniona najprostszym sposobem – wywoływaniem funkcji, która ustawia ten sam seed dla generatora liczb losowych.
setConsistency <- function() {
set.seed(42)
}Ładowanie danych polega na pobraniu ich z internetu, a następnie załadowanie jako tabeli.
loadData <- function(urlPath, localFile = 'life_expectancy.tsv') {
remote_data_file <- download.file(urlPath, destfile = localFile, quiet= TRUE)
return(read.csv(TMP_LOCAL_FILE, quote='"'))
}
FILE_URL <- 'http://www.cs.put.poznan.pl/alabijak/emd/projekt/Life_Expectancy_Data.csv'
TMP_LOCAL_FILE <- 'life_expectancy.csv'
data <- loadData(FILE_URL, TMP_LOCAL_FILE)Podczas wstępnego przetwarzania danych zmieniane są nazwy kolumn, które są niejednolite. Dodatkowo, spożycie alkoholu oraz procent rządowych wydatków na ochronę zostały wyzerowane, jeśli były “NA”, a rok był 2015. Te przykłady nie były usuwane, gdyż każdy przykład nie posiadał tych danych, co może sugerować po prostu na brak danych z tych lat, a większosza ilość danych pomaga klasyfikatorom.
Następnie,przykłady posiadające “NA” w jakimkolwiek atrybucie zostały usunięte, a także utworzony został dodatkowy zbiór posiadający jedynie atrybuty liczbowe.
fixData <- function(df) {
# Rename data
fixedData <- rename(df,
Life.Expectancy = Life.expectancy,
Infant.Deaths = infant.deaths,
Percentage.Expenditure = percentage.expenditure,
Under.Five.Deaths = under.five.deaths,
Total.Expenditure = Total.expenditure,
Thinness.Years.10.19 = thinness..1.19.years,
Thinness.Years.5.9 = thinness.5.9.years,
Income.Resources.Composition = Income.composition.of.resources)
fixedData <- within(fixedData, Alcohol[is.na(Alcohol) & Year == 2015] <- 0)
fixedData <- within(fixedData, Total.Expenditure[is.na(Total.Expenditure) & Year == 2015] <- 0)
return(fixedData)
}
cleanData <- function (df) {
return(na.omit(df))
}
removeChrColumns <- function(df) {
library(dplyr)
return(df %>%
select_if(~!is.character(.)))
}
fixedData <- fixData(data)
cleanedData <- cleanData(fixedData)
numericData <- cleanData(removeChrColumns(cleanedData))Rozmiar danych przedstawiony jest poniżej.
cat(paste('Zbiór danych zawiera', nrow(cleanedData), 'przykładów.\n'))## Zbiór danych zawiera 1777 przykladów.cat(paste('Każdy z przykładow posiada', ncol(cleanedData), 'atrybutów.\n'))## Kazdy z przykladow posiada 22 atrybutów.Do podstawowych statystyk zostało wybrane tylko kilka atrybutów. Jak możemy zaobserwować, średnia długość życia wacha się od 44 do 89 lat, ze średnią na poziomie 69 lat oraz medianą 71 lat. Jest to pozytywna wiadomość, gdyż oznacza, iż krótka długość życia występuje w niewielu krajach. Widać także bardzo duże dysproporcje w liczbie zgonów dzieci oraz na wirusa HIV.
library(dplyr)
library(tidyr)
cleanedData %>%
tibble::as_tibble() %>%
select(Life.Expectancy, Schooling, Infant.Deaths, HIV.AIDS) %>%
summarise(across(everything(), list(mean = mean, median = median, min = min, max = max), .names = "{.col}_{.fn}")) %>%
gather(variable, value) %>%
separate(variable, c("var", "funkcja"), sep = "\\_") %>%
spread(var, value) %>%
relocate(Life.Expectancy, .after = funkcja)## # A tibble: 4 x 5
## funkcja Life.Expectancy HIV.AIDS Infant.Deaths Schooling
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 max 89 50.6 1600 20.7
## 2 mean 69.4 1.90 32.2 12.2
## 3 median 71.7 0.1 3 12.3
## 4 min 44 0.1 0 4.2Przedstawione są histogramy wszystkich atrybutów. Możemy na nich zaobserwować niepokojące ekstrema – wygląda to tak, jakby niektóre z krajów żyły w ciągłej biedzie, a niektóre w ciągłych luksusach (gdyż nie są to rozkłady normalne). W szczególności widać to w GDP, gdzie większość krajów znajduje się bardzo nisko, a tylko elitarnych kilka osiąga duże wartości.
Do porównania korelacji atrybutów został zastosowany współczynnik korelacji Pearsona. Na wykresie poniżej przedstawione są (razem z odbiciem lustrzanym) wszystkie wartości. Warto wspomnieć, iż krzyżyki zaznaczone są tam, gdzie korelacja jest mniejsza niż 0.05 (tz. poziom istotności).
Możemy zauważyć, iż liczba populacji nie wpływa na długość życia. Także rok oraz procent rządowych wydatków nie jest bardzo mało skorelowany z długością życia. Może to być wynikiem różnych krajów “rozwiniętych”, gdzie procent ten jest utrzymywany na stopniu, który zapewnia odpowiedni poziom opieki medycznej (procent może się różnić, ale opieka pozostaje podobna).
library(ggcorrplot)
library(ggplot2)
corr <- round(cor(numericData, method = "pearson"), 1)
p.mat <- cor_pmat(numericData)
dev.new(width=8, height=8)
gg <- ggcorrplot(corr, p.mat = p.mat, hc.order=TRUE,
type='full', tl.cex=8,
sig.level = 0.05, insig= "blank")
fig <- ggplotly(gg)
figŚrednia długość życia dla różnych krajów została przedstawiona na wykresie słupkowym poniżej. Kraj o statusie “Developing” zostały oznaczone na niebiesko, podczas gdy “Developed” – na zielono. Możemy zauważyć, iż kraje rozwinięte są przeważnie w czołówce średniej długości życia, nie znaczy to jednak, iż kraje rozwijające się są gorsze – niektóre z nich, takie jak Kanada i Francja są lepsze od większości krajów rozwiniętych.
library(plotly)
countryNames <- unique(cleanedData$Country)
countryData <- cleanedData %>%
filter(Life.Expectancy != 0) %>%
group_by(Country, Status) %>%
dplyr::summarize(Mean.Life.Expectancy = mean(Life.Expectancy), .groups = 'drop')
cbPalette <- countryData$Status
cbPalette[cbPalette == "Developed"] <- "green"
cbPalette[cbPalette == "Developing"] <- "#56B4E9"
fig <- plot_ly(
countryData,
x = ~Country,
y = ~Mean.Life.Expectancy,
type = 'bar',
marker = list(color = cbPalette)
)
figWytrenowałem regresor – Random Forest Regressor – na wszystkich danych numerycznych tego zbioru. Została zastosowana 2-krotna walidacja krzyżowa z 5-cio krotnym powtórzeniem. Dane zostały przeskalowane.
library(caret)
setConsistency()
inTraining <-
createDataPartition(
y = numericData$Life.Expectancy,
p = .75,
list = FALSE,)
trainingSet <- numericData[inTraining,]
testingSet <- numericData[-inTraining,]
rfGrid <- expand.grid(mtry = seq(2, ncol(numericData) - 1))
gridCtrl <- trainControl(
method = "repeatedcv",
number = 2,
repeats = 5)
fitTune <- train(Life.Expectancy ~ .,
data = trainingSet,
method = "rf",
metric = "RMSE",
preProc = c("center", "scale"),
trControl = gridCtrl,
tuneGrid = rfGrid,
ntree = 10)
fitTune## Random Forest
##
## 1334 samples
## 19 predictor
##
## Pre-processing: centered (19), scaled (19)
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 666, 668, 667, 667, 667, 667, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 2 2.865870 0.8981423 2.008806
## 3 2.611478 0.9138562 1.823620
## 4 2.496997 0.9209252 1.729448
## 5 2.414434 0.9260085 1.651294
## 6 2.334129 0.9300500 1.604995
## 7 2.349975 0.9291291 1.599798
## 8 2.281503 0.9329390 1.568389
## 9 2.303874 0.9315994 1.572533
## 10 2.300910 0.9315208 1.563833
## 11 2.288248 0.9320849 1.549882
## 12 2.274257 0.9330773 1.553570
## 13 2.288128 0.9322114 1.554789
## 14 2.282676 0.9322473 1.545455
## 15 2.297981 0.9313793 1.561153
## 16 2.286173 0.9322201 1.551356
## 17 2.328094 0.9294172 1.578184
## 18 2.300977 0.9314149 1.557778
## 19 2.324170 0.9297220 1.572661
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 12.library(ggplot2)
library(plotly)
gg <- ggplot(fitTune) + theme_bw()
fig <- ggplotly(gg)
figJak widać, regresor zwraca bardzo dobre wyniki, na poziomie 0.95 R2. Oznacza to prawdopodobnie, iż średnia jakość życia może zostać bardzo prosto oszacowana na podstawie danych o kraju.
rfTuneClasses <- predict(fitTune,
newdata = testingSet)
cat("\nTesting set:\n")##
## Testing set:postResample(
as.numeric(rfTuneClasses),
testingSet$Life.Expectancy)## RMSE Rsquared MAE
## 1.9878760 0.9487156 1.2810903Jest to jednakże klasyfikator, które prawdobnie jest bardzo szczegółowy, gdyż ilość wybranych predyktorów wynosi 12. Kolejny regresor jest dokładnie taki sam jak poprzedni, aczkolwiek dane do jego treningu do 6 najbardziej korelujących z długością życia atrybutów.
library(caret)
library(dplyr)
bestCorr <- cor(numericData) %>%
as.data.frame() %>%
dplyr::select(Life.Expectancy) %>%
abs(.) %>%
arrange(desc(Life.Expectancy)) %>%
slice(1:7)
lessData <- numericData %>%
dplyr::select(as.vector(row.names(bestCorr)))
setConsistency()
inTraining <-
createDataPartition(
y = lessData$Life.Expectancy,
p = .75,
list = FALSE,)
trainingSet <- lessData[inTraining,]
testingSet <- lessData[-inTraining,]
rfGrid <- expand.grid(mtry = seq(2, ncol(lessData) - 1))
gridCtrl <- trainControl(
method = "repeatedcv",
number = 2,
repeats = 5)
fitLessTune <- train(Life.Expectancy ~ .,
data = trainingSet,
method = "rf",
metric = "RMSE",
preProc = c("center", "scale"),
trControl = gridCtrl,
tuneGrid = rfGrid,
ntree = 10)
fitLessTune## Random Forest
##
## 1334 samples
## 6 predictor
##
## Pre-processing: centered (6), scaled (6)
## Resampling: Cross-Validated (2 fold, repeated 5 times)
## Summary of sample sizes: 666, 668, 667, 667, 667, 667, ...
## Resampling results across tuning parameters:
##
## mtry RMSE Rsquared MAE
## 2 2.304066 0.9315043 1.587042
## 3 2.263488 0.9334297 1.547365
## 4 2.254369 0.9340318 1.517265
## 5 2.278465 0.9325793 1.527133
## 6 2.291570 0.9318532 1.543398
##
## RMSE was used to select the optimal model using the smallest value.
## The final value used for the model was mtry = 4.rfLessTuneClasses <- predict(fitLessTune,
newdata = testingSet)
cat("\nTesting set:\n")##
## Testing set:postResample(
as.numeric(rfLessTuneClasses),
testingSet$Life.Expectancy)## RMSE Rsquared MAE
## 1.9010933 0.9529699 1.1967020Jak możemy zauważyć, liczba predyktorów zmniejszyła się do 4, a wynik wzrósł. Nie jest to wzrost statystycznie coś mówiący, jednakże sam fakt, iż wynik pozostał na tym samym poziomie 0.95 R2 mówi nam, iż ten klasyfikator powinien być dużo lepszy w generalizacji dla nowych danych.
library(ggplot2)
library(plotly)
gg <- ggplot(fitLessTune) + theme_bw()
fig <- ggplotly(gg)
fig